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INTRODUCTION 


The  relationship  between  free-air  gravity  anomalies  and 
sea-floor  bathymetry  is  one  of  the  r  :>st  important  problems  in 
marine  gravity.  For  purely  practical  reasons,  deciphering 
the  nature  of  such  a  relationship  is  vital  to  predicting 
quantitatively  the  gravity  field  in  oceanic  areas  in  which 
only  bathymetry  is  known.  Closely  related  to  this  practical 
problem  is  the  question  of  the  physical  basis  for  any  observed 
or  proposed  dependence  of  gravity  upon  topography.  Such  a 
question  bears  on  the  mechanisms  for  creation  and  evolution 
of  oceanic  crust  and  lithosphere  and  on  the  posssible  inter¬ 
action  of  the  lithosphere  with  motions  in  the  asthenosphere. 

Marine  gravity  and  sea-floor  topography  cannot  be  related 
by  a  simple  mathematical  expression  (e.g.,  a  linear  relation¬ 
ship)  that  has  validity  in  all  oceanic  environments.  A 
synthesis  of  submarine  gravity  data  and  bathymetry  by  Woolard 
and  Daugherty  (1970)  demonstrated  the  necessity  to  divide 
oceanic  regions  by  tectonic  types  and  the  difficulty  in 
simply  relating  free-air  anomaly  and  depth  even  within  groups 
of  tectonically  similar  environments.  Degree  averages  of 
marine  gravity  anomalies  and  depths  within  selected  regions 
were  also  correlated  with  varying  degrees  of  success  by 
Watts  and  Talwani  (1974) ,  Sc later  et  aT.  (1975) ,  and  Watts 
(1976) . 

A  necessary  step  in  improving  what  heretofore  have 
generally  been  strictly  empirical  attempts  to  derive  a  rule 
relating  gravity  and  bathymetry  over  a  selected  region  is 
to  take  into  proper  account  the  findings  of  plate  tectonic 
studies  of  the  oceans.  At  mid-ocean  ridges,  the  lithosphere 
is  thin,  perhaps  no  thicker  than  the  crust  (Francis  and 
Porter,  1973;  Solomon  and  Julian,  1974;  Orcutt  et  a^.  ,  1975; 
Rosendahl  et^  al.  ,  1976)  ,  and  sea-floor  topography  except 
at  the  shortest  wavelengths  is  isostatical ly  compensated 
with  a  shallow  compensation  depth  (Dorman,  1975;  McKenzie 
and  Bowin,  1976;  Cochran,  1979;  McNutt,  1979).  Sea-floor 
topography  created  by  mid-plate  volcanic  activity  is  com- 

3 


r.iyrwr 


pensated  much  deeper  because  of  the  substantially  thicker 
lithosphere,  and  the  compensation  is  regional  (Vening 
Meinesz,  1941),  involving  flexure  of  the  oceanic  plate 
(Walcott,  1970;  Watts  and  Cochran,  1974;  Watts  et  al. , 

1975,  1980;  Watts,  1978b;  Detridt  and  Watts,  1979) ,  the 
effective  elastic  thickness  of  which  increases  with  sea¬ 
floor  age  (Caldwell  and  Turcotte,  1979).  At  trench  systems, 
topography  is  not  isostatically  compensated  owing  to  the 
large  dynamic  forces  associated  with  lithosphere  subduction 
(e.g.,  Vening  Meinesz,  1954).  Still  controversial  are 
explanations  for  frequent  correlations  in  very  long  wave¬ 
length  gravity  and  topographic  anomalies  in  the  oceans 
(Menard,  1973;  Anderson  et:  aJL.  ,  1973;  Weissel  and  Hayes, 

1974;  Sclater  et  al .  ,  1975;  Marsh  and  Marsh,  1976;  Cochran 
and  Talwani,  1977;  Watts,  1978a) ,  with  asthenospheric  flow 
the  most  exciting  but  still  unproven  hypothesis. 

This  report  gives  the  mathematical  framework  for  a 
simple  relation  between  topography  and  gravity  in  two 
dimensions  for  stable  ocean  basins.  The  conceptual  basis 
for  the  relation  is  flexure  theory  for  thin  elastic  plates 
loaded  from  above.  The  evolution  with  sea-floor  age  of 
lithospheric  temeprature  and  rheology  is  abstracted  to  an 
effective  elastic  layer  thickness  which  grows  with  plate 
age.  In  order  to  derive  gravity  from  topography  using 
this  relation,  both  the  lithospheric  age  and  the  age  of 
any  more  recently  superposed  volcanic  constructs  (islands, 
seamounts,  aseismic  ridges)  must  be  known.  Guidelines  are 
given  for  estimating  these  ages.  The  bathymetry-gravity 
relation  is  tested  using  data  from  the  central  Pacific  basin. 
The  test  is  not  successful,  for  reasons  which  are  described. 
A  number  of  suggestions  are  given  for  future  work. 
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THEORY 


We  seek  a  relation  between  bathymetry  and  gravity  in 
stable  ocean  basins.  We  shall  work  explicitly  in  spatial 
dimensions,  rather  than  use  one-  or  two-dimensional  wave- 
number  representations,  so  as  to  permit  consideration  of 
regions  of  arbitrary  geometry  and  spatial  extent. 

We  shall  assume  that  the  sea- floor  topography  is  com¬ 
posed  of  three  parts:  (1)  the  long-wavelength  deepening 
of  the  sea  floor  with  age  due  to  thermal  contraction  of 
the  lithosphere  (Sclater  ejt  al.  ,  1971);  (2)  the  volcanic 
topography  emplaced  as  a  load  on  top  of  the  lithosphere; 
and  (3)  the  lithospheric  response  to  that  load.  We  assume 
that  the  effect  of  (1)  on  the  bathymetry  may  be  removed  with 
a  suitable  age-depth  relation  (Sclater  <rt  al. ,  1975;  Parsons 
and  Sclater,  1977;  Cochran  and  Talwani,  1977);  its  effect 
on  gravity  is  negligible  far  from  ridges  (Lambeck,  1972) . 

We  discuss  the  validity  of  these  assumptions  further  at  the 
end  of  this  report. 

Let  the  lithosphere  be  modeled  as  a  thin,  spherical, 

elastic  shell  of  thickness  T,  Young's  modulus  E,  and 

Poisson's  ratio  v,  overlying  a  fluid  interior  of  density 

p  .  Let  R  be  the  radius  to  the  midplane  of  the  shall, 
m 

let  g  be  the  gravitational  acceleration  at  radius  r  =  R  -  T/2, 
and  let  D  be  the  flexural  rigidity  of  the  shell 


(1) 


Consider  a  uniform  vertical  load  q  (force  per  unit  area) 
applied  to  a  circular  unit  area  on  the  surface  of  the  shell. 
The  vertical  deflection  (positive  if  downward)  of  the  litho¬ 
sphere  is  given  by 

w  =  - - kei  f;  (2) 

2 7T  (ET/r’  +  Ap  g) 

(Brotchie  and  Sylvester,  1969;  Brotchie,  1971),  where  A  is 


F 


the  distance  from  the  load  center,  normalized  by  the  radius 
of  relative  stiffness 

D _ 

ET/R2  +  Ap  g 

kei  is  a  Bessel-Kelvin  function  of  order  zero  (Abramowitz 
and  Stegun,  1964) ;  and  Ap=  p^  -  p^  where  pw  is  the  density 
of  seawater  if  the  load  is  applied  at  the  sea  bottom  and 
zero  if  the  load  is  applied  on  land.  Note  that  (2)  is 
equivalent  to  the  expression  (Brotchie  and  Sylvester,  1969) 
for  the  response  to  a  point  force  p  since  p  =  q7rd2£2  where 
i  =  d  is  the  radius  of  the  circular  unit  area.  Equations 
(2)  and  (3)  may  be  simplified  because  for  oceanic  litho¬ 
sphere  ET/R2  ~  10 1  and  Ap  g  -  2  x  10 3  in  c.g.s.  units, 

so  Ap  g  >>  ET/R2  and 


w 


q 

2irAp  g 


kei  £. 


(4) 


We  may  generalize  (4)  to  a  distributed  load.  Let 
q (x, y )  be  the  force  per  unit  area  exerted  on  the  litho¬ 
sphere  by  topography.  Then  the  deflection  w  is  given  by 


w (x,y) 


_ 1 _ 

2ff Ap  g 


(x* ,y ' ) 


kei 


( r ' / £ )  dx ' dy ' 


(5) 


where  r'  =  [  (x1  -  x)2  +  (y'  -  y)  2]  .  Equation  (5)  may 

be  thought  of  as  the  double  convolution  of  q  with  the 
flexural  response  function 


<J>  (x,y) 


2irAp  g 


kei 


(x: 


ilL 


(6) 


which  in  turn  is  a  function  of  the  local  flexural  rigidity 
D,  or  equivalently  the  effective  elastic  lithospheric 
thickness  T,  through 


A  single  convolution  of  (two-dimensional)  load  with  response 
function  was  also  used  by  Roufosse  and  Parsons  (1977)  in 
their  study  of  the  Hawaiian  ridge.  If  the  load  is  emplaced 
in  a  time  short  compared  to  the  lithospheric  age  and  if 
viscous  relaxation  of  stress  subsequent  to  initial  litho¬ 
spheric  flexure  has  been  minimal,  then  l  is  simply  a  function 
of  the  time  difference  between  plate  age  and  load  emplacement 
age. 

Unfortunately,  the  load  q(x,y)  is  not  known  a  priori 
but  rather  the  ocean  floor  elevation  h(x,y)  with  respect  to 
some  arbitrary  datum  is  known.  As  noted  above,  h  contains 
contributions  from  both  the  topographic  load  and  the  litho¬ 
spheric  response.  We  may  resolve  this  difficulty,  however, 
by  iteration.  Assume  initially  that  q(x,y)  =  Ap'gh(x,y) 
where  Ap*  =  -  p  and  pfc  is  the  density  of  the  major  units 

constituting  the  topographic  variations.  Then  use  (5)  to 
calculate  w(x,y).  Now  set  q(x,y)  =  Ap'g(h  +  w) .  Recalculate 
w  from  (5) ,  and  repeat  the  last  two  steps  until  q  and  w 
converge  to  a  steady  solution. 

Once  a  self-consistent  decomposition  of  h  into  load 
topography  q/Ap'g  plus  plate  deflection  w  is  achieved, 
the  calculation  of  gravity  is  straightforward.  Two  terms 
contribute  to  the  gravity:  (i)  the  attraction  of  the  topog¬ 
raphy  and  (ii)  the  deflection  of  the  Moho  and  of  any  other 
density  contrast  interfaces  within  the  lithosphere.  Any 
density  contrast  between  the  base  of  the  lithosphere  and 
the  asthenosphere  would  also  contribute  to  (ii)  ,  but  this 
contribution  is  probably  negligible.  For  both  (i)  and  (ii) 
the  gravity  anomaly  may  be  written  in  the  form  of  that  due 
to  a  surface  mass  distribution  o(x,y) 

g(x,y,z)  =  Gz  (  1 - 0  - - dx'dy' 

J  J  [ (x'  -  x) 2  +  (y'  -  y) ?+  z2]3/2 

(7) 

where  G  is  the  gravitational  constant  and  z  >  0  is  the  ver¬ 
tical  distance  between  the  observation  point  (x,y,z)  and  the 


horizontal  plane  on  which  a  is  evaluated.  For  the  topographic 
contribution,  o  =  Ap'  h  on  the  plane  of  the  topographic  datum 
(h  =  0) .  For  the  plate  deflection  contribution,  a  =  Ap  w 
on  the  plane  corresponding  to  mean  Moho  depth. 


IMPLEMENTATION 


The  theoretical  relationship  between  topography  and 
gravity  outlined  above  can  be  applied  to  any  finite  area 
over  which  bathymetry  and  age  are  specified  as  functions 
of  position.  We  have  applied  these  concepts  to  a  region, 
described  in  the  next  section,  in  which  degree  averages 
of  ocean  floor  depth  (corrected  for  mean  age)  and  free  air 
gravity  are  known.  Implementation  of  the  theory  as  a  usable 
algorithm  requires  estimates  of  volcanic  load  ages  and  of 
the  evolution  of  the  flexural  response  function  with  age, 
and  calculational  schemes  for  evaluating  the  convolution 
integrals  in  (5)  and  (7)  . 

Some  care  must  be  exercised  in  deciding  the  age  of  a 
volcanic  construct  relative  to  the  age  of  the  surrounding 
seafloor.  For  islands  this  is  not  a  serious  problem  as 
the  exposed  or  cored  (if  a  coralline  island)  bedrocks  can 
be  dated.  For  seamounts  some  simple  rules  will  be  helpful. 

One  such  rule  is  that  an  inactive  topographic  feature  will 
sink  with  respect  to  sea  level  at  the  rate  of  its  surrounding 
abyssal  sea  floor,  a  rate  which  is  associated  with  thermal 
contraction  of  the  lithosphere  and  which  is  a  well  known 
function  of  ocean-floor  age  (Sclater  et_  al.  ,  1971)  .  On 
this  basis,  for  instance,  the  Ninetyeast  ridge  can  be 
shown  to  have  been  generated  at  the  southeast  Indian  ridge 
(Sclater  and  Fisher,  1974) .  Much  of  the  ocean  floor  topog¬ 
raphy,  in  fact,  was  apparently  generated  at  or  near  ridge 
crests  on  very  young  oceanic  lithosphere  (McKenzie  and  Bowin, 
1976)  . 

The  flexural  rigidity  D  and  the  effective  elastic 

range  thickness  T  are  known  to  vary  from  values  in  the 
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range  3  x  10  -  2  x  10  dyne-cm  and  3-13  km,  respectively, 

near  mid-ocean  ridges  (Cochran,  1979;  Detrick  and  Watts,  1979; 
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McNutt,  1979)  to  values  in  the  range  3  x  10  to  10  dyne-cm 
and  15  -  30  km,  respectively,  for  older  lithosphere  (Walcott, 
1970;  Watts  and  Cochran,  1974;  Watts  et  al.,  1975,  1976; 
Caldwell  et.  al.  ,  1976;  Watts,  1978b;  Suyenaga,  1979).  The 


assumption  that  T  scales  as  an  isotherm  depth  is  consistent 
with  the  data  (Caldwell  and  Turcotte,  1979) ,  so  that  a 
standard  oceanic  plate  thermal  model  (Parsons  and  Sclater, 
1977)  and  a  few  good  measures  of  D  and  T  are  sufficient  to 
define  T  as  a  function  of  lithosphere  age  minus  load  age. 

Evaluation  of  the  necessary  convolution  integrals 


has  been  conducted  as  follows.  For  each  degree  square  we 
first  subtract  the  mean  bathymetric  depth  from  the  depth 
predicted  from  a  spreading  plate  model  for  the  appropriate 
region.  This  residual  depth  is  then  treated  as  a  combin¬ 
ation  of  lithosphere  load  and  response  to  that  load.  For 
the  Green's  function  for  the  loading  problem  we  use  the 
response  of  the  lithosphere  to  a  circular  load  of  radius 


a  such  that 

f 


ira2  equals  the  area  of  the  degree  square: 


-r^—  (a 
Apg 


ker ' 


ber  F  -  ct  kei' 


bei  E,  +  1) 


w  =  < 


for  ^  <  a  =  - 


(8) 


.  —  (ber'  a  ker  £  -  bei'  a  kei  F) 

Apg 


for  £;  >  a  = 


(Brotchie,  1971) ,  where  ker,  ber  and  bei  are  additional 
Bessel-Kelvin  functions  of  zero  order  and  the  prime  denotes 
first  derivative  (Abramowitz  and  Stegun,  1964) .  The 
flexural  length  l  (or  flexural  rigidity  D)  can  be  taken 
appropriate  to  the  lithosphere  at  the  time  of  application 
of  the  load  on  the  degree  square  in  question.  The  convo¬ 


lution  integral  (5)  is  converted  to  a  sum.  Let  w^  be  the 


subsidence  at  the  center  of  the  ith  degree  square,  let  q 


J 


be  the  load  on  the  jth  degree  square  and  let  w^  be  the 
contribution  to  w^  from  q^ .  Then  each  w^  may  be  estimated 
from  (8)  after  setting  F,  to  the  appropriate  (normalized) 
distance  between  the  centers  of  the  ith  and  jth  square,  and 


10 


»• . 


m 


1 


w . 
1 


E 

D 


w .  . 

ID 


(9) 


As  noted  above,  iteration  is  necessary  to  determine  a 

self-consistent  set  of  values  for  q^  and  w^  from  a  given 

set  of  observed  residual  depths.  In  practice,  about  five 

iterations  are  required  before  the  sum  of  the  squared 

differences  in  w.  (or  q.)  between  successive  iterations 
1  9  3 

is  less  than  1  part  in  10  . 

Expression  (7)  for  the  predicted  free-air  gravity 
anomaly  is  evaluated  using  the  Taylor  series  algorithm 
of  Morrison  (1976)  for  calculating  the  gravitational  poten¬ 
tial  from  a  surface  density  distribution  specified  over  a 
latitude-longitude  grid.  The  contribution  from  each  degree 
square  is  of  the  form 


U  (r ,  <t>,  9)  = 


II 


F  d 4> *  d0' 


A<J>'  A0'  [FQ  +  (F( 


A<J>'  2  +  F00  A<T  2)  ]  (9) 


where 


F  = 


G  o  r‘  2  cos 4)' 
R 


<1 


the  observation  point  is  (r ,<P,0)  in  spherical  coordinates 
(radius,  latitude,  longitude),  the  mass  point  is  (r '  ,  <p' ,  O' )  , 
R  is  the  distance  between  the  two  points  and  A<f>’  and  a©’ 
are  the  dimensions  of  the  degree  square.  In  (9),  F0,  F^, 


and  F00  are  F,  62F/6<t>,?  and  <52F/60'2;  respectively,  evaluated 
at  the  center  of  the  di 
obtained  analytically: 


00 

at  the  center  of  the  degree  square.  F  and  F00  may  be 


Gar' 


_  _  2s \n.$-  rr<  [sin4>  cos4>’  -  cos<J>  sin$  '  (cosO  cosG' 

R  R3 

+  sin0  sinO')]  -  -  rr1  [sim)'  sincj>'  +  cos^  cosi) ' 


(cos0  cos6'  +  sin©  sin0 ' )  ]  +  2c_ su  (rr')z  [sin^  cos4> ' 

R5 

-  cost  sin'J>'  (cosO  cosO '  +  sin0  sinO')]2  (10) 
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00 


cos<p' 


rr'  cos<f>  cos <f>*  (cos0  cos0’  +  sin0  sin0') 


Gar' 


(11) 


+  3cos< 


(rr')z  (cosij'  cosi}'1  (-cosO  sinO1  +  sino  cosO1 )  ] 


R2  =  r2  +  r*  2  -  2rr'  [sin4>  sin4>*  +  cos<J>  cos^'  (cos0  cos0'  + 

sin0  sinS1 ) ]  (12) 

Two  separate  contributions  to  U  from  each  degree  block  are 
included  here:  the  contribution  from  topography  (Ap1  = 

1.8  g/cm3  assumed  at  a  datum  plane  5  km  below  sea  level)  and 
the  contribution  from  Moho  deflection  (Ap  =0.6  g/cm3  assumed 
at  a  datum  plane  12  km  below  sea  level) . 

In  practice  the  Taylor  series  truncation  in  (9)  is  a 
good  approximation  only  if  rr'  A^,2/R2  <  1.  This  limits 

(9)  to  elevations  above  the  appropriate  datum  plane  comparable 
to  or  greater  than  the  dimensions  of  the  blocks  used  to  specify 
o .  For  the  degree-square  specification  of  seafloor  depth  and 
free  air  anomaly  in  the  application  considered  here,  we 
calculated  the  potential  at  100  km  altitude  from  (9) ,  and  we 
calculated  the  gravity  anomaly  from  a  centered  finite 
difference  approximation  to  the  radial  derivative  of  the 
potential.  Sea  level  gravity  was  upward  continued  to  100  km 
elevation  for  comparison  with  the  predicted  gravity  anomaly. 

An  alternative  scheme  would  be  to  treat  the  mass  contribution 
(topography  or  Moho  deflection)  degree  square  from  each  as  a 
three-dimensional  prism  and  calculate  the  gravitational 
attraction  using  the  summed  line-integral  method  of 
Talwani  and  Ewing  (1960). 

Fortran  programs  to  calculate  the  load  q  and  subsidence 
w  from  residual  depths  and  to  calculate  the  gravity  from 
the  topography  and  subsidence  are  available  from  the  author. 
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APPLICATION 


We  have  applied  the  bathymetry-gravity  algorithm 
described  above  to  a  portion  of  the  central  Pacific  basin, 
bounded  by  15°  and  35°  N  latitudes  and  150°  and  165°  W 
longitudes.  This  region  was  chosen  because  (1)  gravity, 
bathymetry,  and  lithosphere  age  are  all  well  known;  and 
(2)  there  are  large  superposed  volcanic  loads  on  the  plate 
of  younger  ages,  notably  the  Hawaiian  island  chain.  As  we 
shall  see,  however,  the  selection  of  this  region  in  retro¬ 
spect  was  far  from  ideal  because  of  the  presence  of  sea¬ 
floor  topographic  variations  of  a  type  different  from  those 
assumed  in  the  theory  used  here. 

The  gravity  from  the  region  is  taken  from  the  degree 
averages  of  Watts  and  Leeds  (1977) .  Their  values  for  the 
area  in  question  are  shown  in  Table  1.  The  gravity  contin¬ 
ued  upward  to  100  km  above  sea  level  is  shown  in  Table  2. 

The  bathymetry  is  taken  from  a  1978  NORDA  compilation 
of  degree  averages  provided  by  AFGL  (T.P.  Rooney,  personal 
communication,  1978) .  A  few  blanks  and  obvious  errors  in 
this  data  set  were  filled  in  by  estimating  the  appropriate 
degree  averages  from  charts  of  the  Scripps  Institution 
of  Oceanography  (Chase  et  al. ,  1970) .  The  topographic 
data  are  shown  in  Table  3. 

The  lithosphere  ages  were  taken  from  the  magnetic 
anomaly  maps  of  Pitman  et  jU.  (1974).  Residual  depths  were 
calculated  from  observed  bathymetry  using  the  theoretical 
depth-age  profile  of  Parsons  and  Sclater  (1977)  for  the 
north  Pacific.  Residual  depths  are  shown  in  Table  4. 

The  residual  depths  were  used  to  calculate  the  load  q 
on  the  plate  and  the  resultant  subsidence  w,  following 
the  procedures  outlined  in  the  previous  section.  A  plate 
thickness  T  of  30  km  was  assumed  for  modeling  the  effect 
of  the  Hawaiian  ridge  (Watts,  1978),  and  a  plate  thickness 
of  5  km  was  used  for  modeling  the  effect  of  other  topography 
(assumed  to  have  been  generated  at  or  near  a  spreading 
center)  (Cochran,  1979).  The  resultant  distributions  of 


q  and  w  are  given  in  Tables  5  and  6,  respectively. 

The  gravity  at  100  km  altitude  predicted  from  the 
loading  and  subsidence  model  of  Tables  5  and  6  is  shown 
in  Table  7. 

It  is  immediately  apparent  that  the  predicted  gravity 
does  not  match  the  observations.  The  predicted  gravity 
is  much  too  large,  particularly  in  the  southern  half  of  the 
area  modeled,  and  provides  a  much  poorer  predictor  of  the 
rms  gravity  anomaly  than  does  the  a  priori  assumption  of 
zero  anomaly. 

The  difficulty  lies  with  the  large  depth  anomalies 
shown  in  Table  4.  These  large  positive  (shallow  depth) 
anomalies  over  a  substantial  area  are  treated  as  loads  on 
the  lithosphere  by  the  bathymetry-gravity  algorithm,  whereas 
the  regionally  shallow  seafloor  depth  is  almost  completely 
compensated  by  density  anomalies  most  likely  below  the  litho¬ 
sphere  (Watts,  1976).  The  long-wavelength  gravity  is 
correlated  to  the  long-wave1 ength  residual  depth  anomalies, 
but  the  slope  of  a  linear  fit  between  these  quanitities  in 
the  north  Pacific  is  much  less  (22  mgal/km)  than  in  the 
algorithm  used  here  based  on  a  plate  loading  model.  Thus 
the  failure  of  the  algorithm  for  this  region  of  the  Pacific 
is  due  to  the  large  values  for  residual  depth  that  arise  from 
a  process  not  included  in  the  algorithm,  namely  deep  compen¬ 
sation  of  long-wavelength  topography. 


SUGGESTIONS  FOR  FURTHER  WORK 


In  spite  of  the  wide  disagreement  between  Tables  7  and 
2  (predicted  and  observed  gravity,  respectively),  the 
success  of  the  plate  flexure  model  in  isolated  situations 
over  small  areas  prompts  the  belief  that  this  approach 
toward  relating  gravity  and  bathymetry  has  merit. 

For  regions  similar  to  the  north  Pacific  area  studied 
here,  one  of  two  approaches  should  be  pursued  in  further 
work:  either  (1)  a  relationship  between  gravity  and 

bathymetry  based  on  lithospheric  loading  should  be  sought 
only  for  wavelengths  shorter  than  a  few  hundred  kilometers; 
or  (2)  the  process  that  gives  rise  to  long  (>500  km)  wave¬ 
length  gravity  and  depth  anomalies  should  be  explicitly 
modeled  as  part  of  the  algorithm.  If  approach  (1)  is 
followed,  then  either  small  regions  can  be  treated  in 
isolation  after  'regional'  anomalies  have  been  removed, 
or  the  problem  should  be  conducted  in  the  wavenumber  domain 
and  a  high  pass  filter  applied  to  both  the  gravity  and  depth 
data.  (We  avoided  the  wavenumber  domain  here  to  allow  for 
a  spatially  variable  isostatic  response  function,  but 
perhaps  this  capability  is  a  frequently  unnecessary  luxury). 
Approach  (2)  involves  more  free  parameters  in  the  algorithm 
and  a  much  less  certain  physical  basis  for  modeling;  Watts 
(1976)  has  demonstrated  the  approach  and  the  magnitude  of 
the  effects  for  the  region  studied  here. 
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